Nuclear fusion in dense matter: Reaction rate and carbon burning 
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In this paper we analyze the nuclear fusion rate between equal nuclei for all five different nuclear 
burning regimes in dense matter (two thermonuclear regimes, two pycnonuclear ones, and the inter- 
mediate regime) . The rate is determined by Coulomb barrier penetration in dense environments and 
by the astrophysical S-factor at low energies. We evaluate previous studies of the Coulomb barrier 
problem and propose a simple phenomenological formula for the reaction rate which covers all cases. 
The parameters of this formula can be varied, taking into account current theoretical uncertainties 
in the reaction rate. The results are illustrated for the example of the ^^C+^^C fusion reaction. This 
reaction is very important for the understanding of nuclear burning in evolved stars, in exploding 
white dwarfs producing type la supernovae, and in accreting neutron stars. The S-factor at stellar 
energies depends on a reliable fit and extrapolation of the experimental data. We calculate the 
energy dependence of the S'-factor using a recently developed parameter-free model for the nuclear 
interaction, taking into account the effects of the Pauli nonlocality. For illustration, we analyze the 
efficiency of carbon burning in a wide range of densities and temperatures of stellar matter with the 
emphasis on carbon ignition at densities p > 10^ g cm~^. 

PACS numbers: 25.60.Pj;26.50.+x;97.10.Cv 

I. INTRODUCTION 

We will study nuclear fusion rates of identical nuclei in dense stellar matter. This problem is of utmost importance 
for understanding the structure and evolution of stars of various types. Despite the efforts of many authors the 
theoretical reaction rates are still rather uncertain, especially at high densities. The uncertainties have two aspects. 
The first one is related to nuclear physics and is concerned with the proper treatment of nuclear interaction transitions 
(conveniently described in terms of the astrophysical factor S). The other issue is associated with aspects of plasma 
physics and concerns the proper description of Coulomb barrier penetration in a high density many-body system. We 
will analyze both aspects and illustrate the results taking the carbon fusion reaction as an example. 

Considerable experimental effort has been spent on the study of low energy fusion reactions such as ^^C-f ^^C, to 
investigate the impact on the nucleosynthesis, energy production and time scale of late stellar evolution. Nevertheless, 
it has been difficult to develop a global and reliable reaction formalism to extrapolate the energy dependence of the 
fusion cross section into the stellar energy range. The overall energy dependence of the cross section is determined 
by the Coulomb-barrier tunnel probability. One goal of the present work is to apply the Sao Paulo potential model 
to provide a general description of the stellar fusion processes. This model does not contain any free parameter and 
represents a powerful tool to predict average low energy cross sections for a wide range of fusion reactions, as long 
as the density distribution of the nuclei involved in the reaction can be determined. In this context we also seek to 
introduce a phenomenological formalism for a generalized reaction rate to describe all the regimes of nuclear burning 
in a one-component plasma ion system. In this paper we want to demonstrate the applicability of the method on the 
specific example of ^^C-|-^^C, to evaluate the reliability and uncertainty range of the proposed formalism through the 
comparison with the available low energy data. In a subsequent publication we want to extend the model to multi-ion 
systems with the aim of simulating a broad range of heavy-ion nucleosynthesis scenarios from thermonuclear burning 
in hot stellar plasma, to pycnonuclear burning in high density crystalline stellar matter. 

Carbon burning represents the third phase of stellar evolution for massive stars (Af > 8Mq): it follows helium 
burning that converts He- fuel to ^^C via the triple a process. Carbon burning represents the first stage during stellar 
evolution determined by heavy-ion fusion processes (e.g., Ref. The most important reaction during the carbon 
burning phase is the ^^C-|-^^C fusion 2]; additional processes can be ^^C-|-^^0 and ^^O+^^O, depending on the -'^^C/^^0 
abundance ratio which is determined by the ^^C(q!,7)^^0 reaction rate |3jl4]. The most important reaction branches 
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are ^'^C{^'^C,a)'^°lSle (0=4.617 MeV) and i2c(i2c,p)23Na (0=2.241 MeV). Carbon burning in evolved massive stars 
takes place at typical densities p ~ 10^ g cm"'^ and temperatures T ~ (6 — 8) x 10^ K. 

Carbon burning is also crucial for type la supernovae. These supernova explosions are driven by carbon ignition in 
cores of accreting massive CO white dwarfs The burning process proceeds from the carbon ignition region near the 
center of a white dwarf by detonation or deflagration through the entire white dwarf body. The ignition conditions and 
time scale are defined by the ^^C+^^C reaction rate, typically, at T ^ (1.5 — 7) x 10^ K and p ^ (2 — 5) x 10^ g cm^^ 
Depending on the ^^O abundance, other fusion reactions may also contribute. For these high densities the 
reaction cross sections are affected by strong plasma screening, which reduces the repulsive Coulomb barrier between 
interacting ^^C or ^^O nuclei (e.g., Refs. ^8, 9]; also see Section UTTj) . 

Explosive carbon burning in the crust of accreting neutron stars has recently been proposed as a possible trigger and 
energy source for superbursts 0,0,^3- ^^^^ scenario small amounts of carbon (3%-10%), which have survived in 
the preceding rp-process phase during the thermonuclear runaway, ignite after the rp-process ashes are compressed by 
accretion to a density p ~ 1.3 x 10^ g cm~^. The ignition of a carbon flash requires an initial temperature T > 10^ K, 
triggering a photodisintegration runaway of the rp-process ashes after a critical temperature T ^ 2 x 10^ K is reached 
[iSj. For these scenarios, carbon burning proceeds in the thermonuclear regime with strong plasma screening (see 
Section ITnl for details). 

At high densities and/or low temperatures the thermonuclear reaction rate formalism is insufficient since the fusion 
process is mainly driven by the high density conditions in stellar matter f Sections IIIII and IIVII - This is particularly 
important for nuclear fusion in the deeper layers of the crust of an accreting neutron star |l4| . At sufficiently high 
p and low T nuclei form a crystalline lattice. Neighboring nuclei may penetrate the Coulomb barrier and fuse owing 
to zero-point vibrations in their lattice sites. In this pycnonuclear burning regime the reaction rate depends mainly 
on the density and is nearly independent of temperature (e.g. 15, 16] ). Pycnonuclear burning regimes may not be 
limited to carbon induced fusion reactions only, but may be driven by a broad range of fusion reactions between stable 
and neutron rich isotopes 'l4l. 

In the following Section^ we discuss the theory of fusion cross sections and calculate the astrophysical 5-factor for 
carbon burning in the framework of a generalized parameter-free potential model. In Sect ion Hill we study the Coulomb 
barrier problem for identical nuclei, and propose an expression for the reaction rate which describes all the regimes 
of nuclear burning in a dense one-component plasma of atomic nuclei. In Section llVI we analyze, for illustration, the 
main features of ^^C burning from high temperature gaseous or liquid plasma to high density crystalline matter. We 
summarize and conclude in Section Ivl 



II. FUSION CROSS SECTION AND ASTROPHYSICAL 5-FACTOR 



Nuclear reactions are possible after colliding nuclei tunnel through the Coulomb barrier. Recently, a parameter- 
free model for the real part of the nuclear interaction (Sao Paulo potential) based on nonlocal quantum effects was 
developed [H El IH • In previous work 21], this model was applied to the study of fusion processes using the 
barrier penetration (BP) formalism for about 2500 cross section data, corresponding to approximately 165 different 
systems. Within the nonlocal model, the bare interaction Vjv (?',£') is connected with the folding potential Vpir), 

VN{r,E)^VF{r)e-^-"l-\ (1) 

where c is the speed of light, E is the particle collision energy (in the center-of-mass reference frame), v is the local 
relative velocity of the two nuclei 1 and 2, 

v2 (r, E)^-[E-Vc{r)-VN{r,E)] , (2) 

Vc{r) is the Coulomb potential, p = AiA2mu/{Ai + A2) is the reduced mass, and is the atomic mass unit. The 
folding potential depends on the matter densities of the nuclei involved in the collision: 

Vf{R) = J Pi(ri) p2(r-2) Vq d{R - n + r^) dr^ (3) 

with Vb = —456 MeV fm"^. The use of the matter densities and delta function in Eq. Q corresponds to the zero- 
range approach for the folding potential, which is equivalent [23 to the more usual procedure of using the M3Y 
effective nucleon-nucleon interaction with the nucleon densities of the nuclei. The advantage in adopting the Sao 
Paulo potential to describe the fusion cross section relies on the fact that no additional parameter is necessary once 
the density distribution of the participating nuclei has been determined. The model is therefore a good choice for a 
generalized treatment of low energy heavy- ion fusion reactions. 
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There are several ways to determine the nuclear density distribution [20|, |23, |2J, |25|, |27j, |28|. Density 
Functional Theories (DFT) provide for example a successful description of many nuclear ground state properties, 
in particular, of charge distributions in the experimentally known region. Since these theories are universal in the 
sense that their parameter sets are carefully adjusted and valid all over the periodic table, one can expect that they 
also yield reliable predictions for nuclei far from stability. Non-relativistic density functionals, such as the Skyrme or 
Gogny functional, have been widely used in the literature. In recent years, relativistic density functionals have played 
an increasingly important role since they provide a fully consistent description of the spin-orbit splitting. This is of 
greatest importance for nuclei far from stability. The spin-orbit splitting determines the shell structure, the most 
basic ingredient in any microscopic theory of finite nuclei. In fact, the results obtained with relativistic functionals are 
in very good agreement with experimental data, throughout the periodic table, despite having a smaller number of 
adjustable parameters in comparison with the non-relativistic case. Best known is the relativistic Hartrcc-Bogoliubov 
theory J23, ^2^, '24] , which includes pairing correlations with finite range pairing forces. It provides a unified description 
of mean- field and pairing correlations in nuclei. 

These functionals contain a strong density dependence, either through non-linear coupling terms between the 
meson fields (e.g., in the Lagrangians with the parameter sets NL2 ^23\ and NL3 '2^), or by using an explicit density 
dependence for the meson-nucleon vertices (e.g., in the parameter sets DD-MEl 27] and DD-ME2 [28|). 

In the present paper, we consider only spherical nuclear shapes. Pairing correlations are in principle included, but 
they vanish for the ^^C nucleus. In Fig. ^ we compare the calculated densities with experimental data 29]. The 
RHB calculations arc in good agreement with surface properties best described by the NL2, DD-MEl, and DD-ME2 
effective interactions. 

To apply the BP model for calculating fusion cross sections one needs the effective potential defined as a sum of 
the Coulomb, nuclear and centrifugal components: 

Kff ir,E) = Vcir) + VN{r,E) + ^J^±l^ . (4) 

Following the BP model one can associate the fusion cross section with the particle flux transmitted through the 
barrier. 

It is important to point out that the sum in Eq. (jSJ is performed up to a maximum £ wave (icr), which corresponds 
to the greatest value of angular momentum that produces a pocket (and a barrier) in the corresponding effective 
potential, Eq. For waves with effective barrier heights Vse < E, the shape of the effective potential can be 
approximated by a parabola with curvature defined as 



/i dr^ 



(6) 



where Rbi is the barrier radius. In such cases, the transmission coefficients have been obtained through the Hill- 
Wheeler formula 1301: 



= n + exp 



2n{VBi - E) 



(7) 



On the other hand, for £-waves with Vse > E, instead of the Hill- Wheeler formula, we employ a more appropriate 
heuristic treatment based on a WKB approximation l3lj ]: 

T, = [l + exp(^,)]-\ (8) 



where ri and r2 are the classical turning points. At low energies, the WKB method gives values for the transmission 
coefficients which are quite different from those of the Hill- Wheeler formula. In this case, we define the barrier 
curvature by connecting Eqs. (6) and (7): 

2TTiVBi - E) 

'St ■ 
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The overall results provided by the BP model are in very good agreement with the fusion data for energies above 
the s-wave barrier height. For light systems {fi < Sm^) the model also shows very good agreement with fusion data 
at sub-barrier energies |23|- Therefore, the use of the BP model in calculating the fusion cross section at energies of 
astrophysical interest for the ^^C+^^C system is entirely justified. 

Historically, reaction cross sections cr{E) at very low energies, typical for astrophysical conditions, have been 
expressed in terms of the astrophysical S*- factor (e.g., Ref. ,32]), 

S'(£;) = fj (£:)£: e^'^" , (ii) 



where 77 = {ZiZ2e^ /h) ^fi/{2E) is the usual Gamow parameter. 

Considerable efforts have been made over the last decades to measure the ^^C+^^C fusion cross section at very 
low energies [23, 0, IsS EE 1^ Ell • The experimentally determined 5- factors are shown in Fig. El For reaction 
rate calculations the experimental S'-factor needs to be extrapolated towards the stellar energy range, the Gamow 
window, which depends sensitively on the temperature and density conditions of the stellar environment. The typical 
range of energy E for thermonuclear carbon burning, in the center-of-mass reference, varies from 1 to 4 MeV. For 
pycnonuclear carbon burning in the neutron star crusts the energies can be as low as 10 keV. Large discrepancies 
between the different experimental results at low energies complicate a reliable extrapolation of S{E) towards such 
low E. In addition, the S'-factor shows pronounced resonant structures, presumably resulting from quasimolecular 
doorway states. Theoretical calculations of S{E) using the effective interactions NL2, NL3, DD-MEl, DD-ME2 agree 
reasonably well, within a factor ^ 3.5 in the limit E (Fig. [2). Furthermore, the resonant behavior of the 
data cannot be described with the BP calculations because the effects of nuclear structure were neglected. However, 
an average description of the data (neglecting resonant oscillations) for the sub-barrier region {E < 6.0 MeV) is 
reproduced satisfactorily. Such a description of an average S'-factor is quite sufficient since the reaction rate formalism 
relies on the average S-factor behavior over the entire Gamow range. 

In this context it is important to emphasize that the main purpose of this paper is not to investigate the oscillations 
in the ^^C-|-"'^^C fusion excitation function. In order to reproduce the resonances we could, for example, use the concept 
of internal and barrier waves based on a semiclassical description |39l | or adopt the R- matrix formalism (e.g., Michaud 
and Vogt ^3)- However, neither theoretical approach would allow us to extrapolate with confidence the fusion cross 
section to the energy region of astrophysical interest. 

In order to calculate the carbon burning rate, we use the values of S{E) obtained on the basis of the well established 
NL2 effective interaction. As one can see from Fig. the ^^C density distribution obtained using the parameter set 
NL2 can describe satisfactorily the surface properties, which is the most important region for the fusion process at 
low energies. The values of S{E) calculated at i? < 19.8 MeV can be fitted by an analytic expression 

C o 7^0.308 ^ 

S(E) = 5.15 X 10^^ exp <^ -0.428 E ^, } MeV barn, (12) 

where the center-of-mass energy E is expressed in MeV. The formal maximum fit error, 16%, occurs a,t E — 5.8 MeV. 
However, let us bear in mind that the values of S{E) provided by the NL2 model and given by Eq. (|12|l are actually 
uncertain within a factor of 3.5. 

With the aim of investigating the validity of our assumption for the real part of the nuclear interaction, we performed 
an optical model (OM) analysis of the ^^C-|-^^C elastic scattering data at energies around and slightly above the 
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Coulomb barrier |41| . We defined the imaginary part of the optical potential, which accounts for the nuclear absorption 
process, as 

W{r,E) = N,VN{r,E), (13) 

where VN{r,E) is described by Eq. (Q) and Ni = 0.78 was determined by adjusting thirty elastic scattering angular 
distributions corresponding to seven different heavy- ion systems and measured in a very wide energy range |42j| . 
Figure O illustrates a comparison between our OM analysis and five elastic scattering angular distribution data of the 
12(-<_|_12q system. As one can note, it is possible to obtain a reasonable description of the data by adopting the Sao 
Paulo potential to account for the real part of the nuclear interaction, combined with a simple model to describe the 
imaginary part of the optical potential. This means that both elastic scattering and fusion processes can be described 
by the same real part of the nuclear interaction, which has been well accounted by the Sao Paulo potential, Eq. JQ. 
As discussed in Ref. |42j , details on the absorption part of the interaction are not very important for describing the 
elastic scattering data, which allows us to get reasonable estimates for the ^^C-|-"'^^C system. 

Further experiments at lower energies are necessary to confirm the validity of the predicted ^^C+^^C S'-factor and 
its impact on the reaction rate. However, the S-factor is not the only uncertainty for a reliable description of the 
^^C-|-^^C fusion process in stellar matter. The reaction rates are also uncertain because of the problems in calculating 
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the probability of Coulomb barrier penetration in a dense many-body environment. We shall discuss these problems 
in Section 11111 and show that the associated uncertainties are higher than the current uncertainties in the values of 
SiE). 

III. NUCLEAR FUSION RATE IN DENSE MATTER 
A. Physical conditions and reaction regimes 

In the following we will turn our attention to the plasma physics aspects of nuclear burning in dense matter. We will 
focus on the formalism of fusion reactions between identical nuclei {A, Z) + {A, Z) in the wide domain of temperatures 
T and densities p, characteristic for the range of stellar environments outlined above. 

As an example, we consider the ^^C+^^C reaction in stellar matter at conditions displayed in the p—T phase diagram 
in Fig.^ Under these conditions, carbon is fully ionized (either by electron pressure and/or by high temperature) and 
immersed in an almost uniform electron background. The electrons are typically strongly degenerate; their degeneracy 
temperature Tp is shown in the figure. 

The state of ions (nuclei) is determined by the Coulomb coupling parameter F = Z^e^/ (aT), where a — [3/(47rrii)]^/'^ 
is the ion-sphere radius and Ui is the number density of ions; the Boltzmann constant is set fee = 1- If T < 1 (which 
happens at T > T; = Z^e^/a, see Fig. 0J, the ions constitute a Boltzmann gas, while at higher F they constitute 
a strongly coupled Coulomb liquid. The gas transforms smoothly into the liquid, without any phase transition. 
At small T (large F) the liquid can solidify. In the density range displayed in Fig. ^ the solidification occurs at 
T = Tm — ^^e^/aF,n, where F,n = 175 (e.g., De Witt et al. (43|)- The important measure of quantum effects in 
ion motion is provided by the ion plasma frequency ujp — y^4TrZ^e^ni/m or the associated ion plasma temperature 
Tp — hujp (to being the ion mass). As a rule, the quantum effects are strongly pronounced at T below Tp. 

Figure 01 shows that the ion system can have very different properties, depending on T and p. As a result, there 
are five qualitatively different regimes of nuclear burning in dense matter (Salpeter and Van Horn [H). These are 
(1) the classical thermonuclear regime; (2) the thermonuclear regime with strong plasma screening; (3) the thermo- 
pycnonuclear regime; (4) the thermally enhanced pycnonuclear regime; and (5) the zero-temperature pycnonuclear 
regime. The regimes differ mainly in the character of the Coulomb barrier penetration of reacting nuclei. The 
penetration can be greatly complicated by Coulomb fields of ions which surround the reacting nuclei. These fields are 
fluctuating and random (e.g., Alastuey and Jancovici 44] ). 

A strict solution of the barrier penetration problem should imply the calculation of the tunneling probability in 
a random potential, with subsequent averaging over an ensemble of random potentials. This program has not been 
fully realized so far. The exact theory should take into account a range of effects which can be subdivided (somewhat 
conventionally) into classical and quantum ones. The classical effects are associated with classical motion of plasma 
ions and with related structure of Coulomb plasma fields (including spatial and temporal variability of these fields). 
The quantum effects manifest themselves in ion motion (e.g., zero-point ion vibrations), quantum "widths" of ion 
trajectories during Coulomb barrier penetration, and quantum statistics of reacting nuclei. The effects of quantum 
statistics are usually small due to the obvious reason that quantum tunneling lengths are typically much larger than 
nuclear radii. The smallness of these effects has been confirmed by Ogata |43 in path-integral Monte Carlo (PIMC) 
simulations. 

The reaction rates in the classical thermonuclear regime are well known (e.g.. Fowler, Caughlan, and Zimmerman 
[3^'): they have been tested very successfully by the theory and observations of the evolution of normal stars. This 
theory will be only shortly reviewed in the following section. The reaction rates in other regimes have been calculated 
by a number of authors in different approximations. In the following we summarize the main results published after 
the seminal paper by Salpeter and Van Horn '151 (^^^ that paper for references to earlier works). Let us stress that the 
reaction rate is a rapidly varying function of plasma parameters. In the most important density-temperature domain 
it varies over tens orders of magnitude fSection lIV|l . In this situation, a very precise calculation of the reaction rate 
is very difficult but not required for many applications. 

B. Classical thermonuclear reaction rate 

The classical thermonuclear regime takes place at sufficiently high T and low p so that the ions constitute a 
Boltzmann gas (T ^ T;, Fig. 2)). The tunnel probability (penetrability) through the Coulomb barrier depends on 
the energy of the interacting ions; the main contribution to the reaction rate comes from ion collisions with energies 
approximately equal to the Gamow peak energy i?pk (that is much higher than T). This regime is typical for all 
nuclear burning stages in "normal" stars (from the main sequence to pre-supernovae). 
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The thermonuclear reaction rate is expressed by 



where i?pk — Tt/S is the Gamow peak energy and 

2Tn- ) =[ ATH- ) ^''^ 
is the parameter which characterizes the penetrabihty exp(— r). The parameter r can be rewritten as 

r = 3 (^/2)2/3(i?jT)i/3, = mZ^eyff. (16) 
Now the reaction rate can be presented as 

i?th = ^ S{Ep^) Pth Fth, (17) 

where h/{mZ'^e'^) is a convenient dimensional factor, i<th is the exponential function, and Pth is the pre-exponent: 

87ri/3 /E \^^^ 

F.,=exp(-.), . (18) 

The classical thermonuclear reaction rate decreases exponentially with decreasing T. 

C. Thermonuclear burning with strong plasma screening 

The thermonuclear regime with strong plasma screening operates in a colder and denser plasma {Tp < T < T;), 
where ions constitute a strongly coupled classical Coulomb system (liquid or solid). The majority of ions in such a 
system are confined in deep Coulomb potential wells (Z^e^/a > T). The main contribution into the reaction rate 
comes from a small amount of higher-energy, unbound ions with E « iJpk S> T from the tail of the Boltzmann 
distribution. The plasma screening effects are produced by surrounding plasma ions and simplify close approaches of 
the reacting nuclei, required for a successful Coulomb tunneling. This enhances the reaction rate with respect to that 
given by Eqs. (O and lfTH|) . 

The enhancement has been studied by a number of authors, beginning with Salpeter Q; it can reach many orders 
of magnitude. Calculations show that the equations given in Section IIII Bl remain valid in this regime, but the 
penetrability function Fth has to be corrected for the screening effects: 

Fth = F,, ■ exp(-r) F,, = exp{h), (19) 

where Fgc is the enhancement factor, and /i is a function of plasma parameters. 

Plasma screening effects are usually modeled by introducing a mean- force plasma potential H{r). In this approx- 
imation, the reacting nuclei move in a potential W{r) = Z'^e'^/r ~ H{r). The mean-force plasma potential H{r) is 
static and spherically symmetric. It cannot take into account dynamical variations of plasma microfields and their 
instantaneous spatial structures in the course of an individual tunneling event. In the mean- force approximation, the 
function h consists of two parts, h = Hq + hi, where the leading term ho = H{0)/T (^ is calculated assuming 
a constant plasma potential H{r) = H{0) during the quantum tunneling, while hi is a correction due to a weak 
variation of H{r) along the tunneling path. Note that according to simple estimates (e.g., Ref. ?4f3j) typical timneling 
lengths of reacting ions in the thermonuclear regime (where T > Tp) are considerably smaller than the ion sphere 
radius a, and typical tunneling times are much smaller than the plasma oscillation period ^ ^p^- This justifies the 
assumption of almost constant and static plasma potential during a tunneling event. 

The mean-force plasma potential H{r) for a classical strongly coupled system of ions (liquid or solid) can be 
determined using classical Monte Carlo (MC) sampling (e.g., DeWitt et al. 01 )• sampling gives the static 

radial-pair distribution function of ions g{r) = exp{—W{r)/T) which enables one to find H{r). In this way one can 
accurately determine g{r) and H{r) at not too small r (typically, at r > a), because of poor MC statistics of close 
ion separations. The potential H{r) at small r, required for a tunneling problem, is obtained by extrapolating MC 
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values of H{r) to r — > 0; the extrapolation procedure is a delicate subject and may be ambiguous (as discussed, e.g., 
by Rosenfeld (4^1). 

It is only H{0) which is required for finding Hq. For a classical ion system, H{Q) can be determined by H{Q) = AjF, 
where AjF is a difference of Coulomb free energies (for a given system and for a system with two nuclei merging into 
one compound nucleus; e.g., DeWitt et al. (47i|). In this approximation, the enhancement factor of the nuclear reaction 
becomes a thermodynamic quantity and acquires a Boltzmann form, exp(ft,o) = exp(AjF/T), showing that plasma 
screening increases the probability of close separations (and subsequent quantum tunneling) ; Hq becomes the function 
of one argument T. Assuming a linear mixing rule in a multi-component strongly coupled ion system, Jancovici |49| 
obtained ho = 2/o(r) — /o(2^/'^r), where /o(r) is a Coulomb free energy of one ion in a one-component plasma of 
ions (in units of T). In a Coulomb liquid at F > 1 the linear mixing rule is highly accurate (DeWitt and Slattery 
[S^l); the function /o(F) is now determined from MC samphng with ver y h igh accuracy (e.g., Ref. [13,1^). In this 
way the function ho{T) has been calculated in many papers (e.g.. Kefs, lillliliiil), and the results are m very 
good agreement. Let us present the analytical approximation of /io(r) which follows from the recent MC results of 
DeWitt and Slattery 0] for a Coulomb liquid at 1 < F < 170: 

ho = 1.0563 F -I- 1.0208 F"-323i_Q 274g In F- 1.0843. (20) 

However, this accurate expression is inconvenient for further use, and we propose another fit 

ho = C^c T^^V[{C,c/V3)^ + (21) 

where Csc = 1.0754. It approximates e''" with the maximum error of ~40% at F = 170, quite sufhcient for our 
purpose. There may be still some uncertainty of the reaction rate associated with the choice of Csc but it seems to be 
not higher than the uncertainty in the 5- factor (Section^. Our fit function in Eq. H21() is chosen in such a way to 
reproduce also the well known expression ho — > •\/3F'^/^ derived by Salpeter 8] for the classical thermonuclear regime 
(F ^ 1), where ho 1 and the plasma screening is weak. In the Coulomb liquid, at F > 1, we have actually the 
linear function ho = CgcF. Such a function was obtained by Salpeter using a simple model of ion spheres (with 
shghtly lower coefficient, Cg^^P = 1.057). 

Some authors calculated ho and the associated enhancement factor e'*" by extrapolating MC H{r) to r ^ (as 
discussed above). In particular, Ogata et al. 0| employed this formalism to study the enhancement of nuclear 
reactions in one-component and two-component strongly coupled ion liquids. The enhancement factor e''" for a one- 
component ion liquid, calculated in these papers (e.g., Eq. (6) in Ref. 53]), is systematically higher than the factor 
given by Eq. (|20|l or H21() . The difference reaches a factor of approximately 40 for F ^ 170. Because the enhancement 
factor itself becomes as high as e'*" ~ 10^^ at F ~ 170, such a difference is insignificant for many applications. As 
shown by Rosenfeld jS^, the difference comes from the problems of extrapolation of H{r) to r ^ in Refs. [sM l53|. 
The function ho was also calculated by Ogata using direct PIMC method. His result (his Eq. (19)) is in much 
better agreement with Eq. H20I) . The maximum difference of e'^" reaches only a factor of approximately 6 at F = 170. 
Recently new PIMC calculations have been performed by Pollock and Militzer [5^ but the authors have not calculated 
directly /io(F). 

Let us emphasize that the enhancement factor e''"^ derived in a constant mean- force plasma potential H{0), is 
invariant with respect to the order of the mean-force averaging and the tunneling probability calculation. One can 
consider a real (random) plasma potential, constant over a tunneling path in an individual tunneling event. Calculating 
the tunneling probability and averaging over an ensemble of realizations of plasma potentials, one comes (e.g., Ref. 
to the same expression for ho as given by the mean- force potential. 

In addition to e''", the enhancement factor Fgc in Eq. (|19|l contains a smaller factor e'*^, associated with variations of 
the plasma potential along the tunneling path. Numerous calculations of hi have commonly employed the mean-force 
potential H{r). The results are sensitive to the behavior of H{r) at small r (where this behavior is not very certain). 
For example, Jancovici 49] got hi = -(5/32)F(3F/t)2. Note that for the thermonuclear burning (T > Tp, Section 
Irnlll . the ratio 3F/t « n/a - {T/Tpf/^ can be regarded as a small parameter (rt being the tunneling length). It 
is possible that the mean-force approximation is too crude for calculating hi. For that reason, we will not specify hi 
in this section. Our final expression for the reaction rate will include hi, but phenomenologically, when we combine 
reaction rates in all regimes fSection IlII G|l . 

D. Zero-temperature pycnonuclear fusion 

The zero-temperature pycnonuclear regime operates in a cold and dense matter (T well below Tp) in a strongly 
coupled quantum system of nuclei. In this regime the Coulomb barrier is penetrated owing to zero-point vibrations of 
neighboring nuclei which occupy their ground states in a strongly coupled system. One usually considers pycnonuclear 
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reactions in a crystalline lattice of nuclei but they are also possible in a quantum liquid. The main contribution into the 
reaction rate comes from pairs of nuclei which are most closely spaced. The reaction rate is temperature-independent 
but increases exponentially with increasing density as we will discuss in the following. 

Pycnonuclear reaction rates between identical nuclei in crystalline lattice have been calculated by many authors 
using different approximations. In analogy with Eq. (|17() , the resulting reaction rates can be written as 

nf h 

^pyc — ~^ S{Epi^) ^^^2g2 ^py^ -^pycj (22) 
where i^pyc and i^pyc depend on the density and have the form 

i^pyC = exp (-Ccxp/Va) , Ppyc = 8 Cpyc U-blb/X^^K (23) 

The dimensionless parameters Coxp, Cpi and Cpyc are model dependent (see below). The dimensionless parameter A 
is expressed in terms of the mass fraction Xi contained in atomic nuclei (in a one-component ion plasma under study) 
and the mass density p of the medium 

/n^x 1/3 ^J_(l_ pX, 
~ raZ'^e^ V 2 / ~ AZ"^ \A 1.3574 x 10" g cm-3 

For densities p lower than the neutron drip density (^^ 4 x 10"'^^ g cm~^; e.g., Ref. H^), one can set Xi = 1, while for 
higher p one has Xi < I because of the presence of free (dripped) neutrons. 
The reaction rate can be expressed numerically as 

i?pyc = pX, AZ^ S{Epy) Cpyc 10^6 (^_Cexp/\/A) s"! cm-^, (25) 

where p is in g cm~'^ and S(Epw) is in MeV barn. The typical energy of the interacting nuclei is Ep\^ ^ fkOp. 

Table|2lists the values of Cexp, Cpi and Cpyc reported in the literature for two models (1 and 2) of Coulomb barrier 
penetration by Salpeter and Van Horn for six models (3-8) by Schramm and Koonin 0, and for one model (9) 
by Ogata, lyetomi and Ichimaru |53j | . The corresponding carbon burning rates are plotted as a function of density in 
Fig. El In this figure (as well as in Figs. 0] and O we use the astrophysical factors given by the fit expression l|12|l . 
Actually, the ^'-factors are uncertain within one order of magnitude (Section but we ignore these uncertainties 
(because they seem to be much lower than those associated with the Coulomb barrier problem). 

All the authors cited above have treated quantum tunneling by fixing the center-of-mass of reacting nuclei in its 
equilibrium position. All models, except for models 5-8, focus on nuclear reactions in the body-centered cubic (bcc) 
lattice of atomic nuclei. This lattice is thought to be preferable over other lattices, particularly, over the face-centered 
cubic (fee) lattice. The main reason is that the bcc lattice is more tightly bound in the approximation of a rigid 
electron background. However, the difference in binding energies of bcc and fee lattices is small (see, e.g., Ref. 16J), 
and a finite polarizability of the electron background complicates the problem 56] . Therefore, one cannot exclude 
that the lattice type is fee. 

Salpeter and Van Horn calculated the quantum tunneling probability of interacting nuclei in a bcc Coulomb 
lattice using the three-dimensional WKB approximation (most adequate for the given problem) . The authors employed 
two models, static and relaxed lattice (models 1 and 2 in Table P), to account for the lattice response to the motion 
of tunneling nuclei. The static lattice model assumes that surrounding nuclei remain in their original lattice sites 
during the tunneling process. The relaxed lattice model assumes that the surrounding nuclei are promptly rearranged 
into new equilibrium positions in response to the motion of the reacting nuclei. Simple estimates show that the 
actual tunneling is dynamical (neither static not relaxed). Thus, the static-lattice and relaxed-lattice models impose 
constraints on the actual reaction rate. In Ref. |0| the screening potential for the relaxed-lattice model was calculated 
approximately; the energy difference between the initial and fused states was evaluated by subtracting the energies 
of the corresponding Wigner-Seitz (WS) spheres. The relaxed lattice simplifies Coulomb tunneling and increases the 
reaction rate with respect to the static lattice (ef. curves 1 and 2 in Fig.jS)). 

Schramm and Koonin applied this treatment to the bcc and fee, static and relaxed lattices in the same WKB 
approximation. They calculated the screening potential for the relaxed-lattice model with improved accuracy (model 
3 of Table m for the bcc lattice and model 7 for fee). For comparison, they also used the screening potential for 
the relaxed-lattice obtained in the WS approximation (as in Ref. ^3). Unfortunately, they calculated the tunneling 
probability neglecting the correction for to the "curvature of trajectories" of reacting ions. This is the main reason 
for the formal disagreement between the results of Salpeter and Van Horn '15] and Schramm and Koonin 16] for 
the static-lattice and relaxed-lattice- WS models (1 and 2) of bcc crystals. The inclusion of the curvature correction 
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should reduce the constant Cpyc in Eq. (|25|l and the reaction rates calculated in Ref . |T6| . Fortunately, this correction 
can be extracted by comparing Eq. (38) of Ref. with Eq. (31) of Ref. jl^J. In this way we get e^ — 0.067 for the 
static bee lattice, and = 0.050 for the relaxed- WS bcc lattice. After introducing this correction into the coefficients 
Cpyc, obtained formally from the results of Schramm and Koonin, these coefficients become identical to those given by 
Salpeter and Van Horn. Thus, Schramm and Koonin actually exactly reproduce models 1 and 2 of Salpeter and Van 
Horn. The curvature correction for the models 3-8 of Schramm and Koonin have not been determined. We expect it 
to be ~ 0.050 for model 3 and e^ = 0.067 for model 4 (bcc crystals), and we introduced such corrections in Table 
U We introduced, somewhat arbitrarily, the correction e^ = 0.05 in all fee lattice models 5-8. 

The two versions of the screening potential for the relaxed lattice (WS and more accurate) almost coincide. Ac- 
cordingly, models 2 and 3 yield almost the same reaction rates for the bcc lattice, while models 6 and 7 yield nearly 
identical rates for the fee lattice (Fig.O. 

Schramm and Koonin [l^ also took into account the dynamical effect of motion of the surrounding ions in response 
to the motion of tunneling nuclei in the relaxed lattice (models 4 and 8). This effect was described by introducing the 
effective mass of the reacting nuclei. The effective mass appears to be noticeably higher than the real nucleus mass, 
reducing the tunneling probability. It turns out that the reduction almost exactly compensates the increase of the 
tunneling probability due to the lattice relaxation neglecting the effective mass effects jl^]. Accordingly, model 4 of 
Schramm and Koonin [l^ gives almost the same reaction rate as model 1 (for bcc) ; and model 8 gives almost the same 
rate as model 5 (for fee). This means that the two limiting approximations, the static-lattice and relaxed-lattice, yield 
very similar reaction rates. It is natural to expect that the actual reaction rate (to be calculated for the dynamically 
responding lattice) would be the same, and the problem of dynamical tunneling is thus solved 0| . Notice that this 
conclusion is made using the curvature corrections adopted above (whereas the accurate curvature correction for 
the effective mass model has not been calculated). 

Zero-temperature pycnonuclear reactions in bcc crystals were also studied by Ogata, lyetomi, and Ichimaru jH^ and 
Ichimaru, Ogata, and Van Horn j53| using MC lattice screening potentials. These authors considered one-component 
and two-component ion systems. Model 9 of Table represents their results for one-component bcc crystals. In order 
to calculate the tunneling probability, the authors used the mean-force plasma screening potential H{r), obtained from 
MC sampling in a classical bcc crystal at r > a and extrapolated to r ^ (see Section llH C|) . This potential is static 
and spherically symmetric. It cannot take into account the dynamics of lattice response and the anisotropic character of 
the real screening potential in a lattice. Furthermore, the barrier penetration was calculated by solving numerically the 
effective radial Schrodinger equation. This procedure is more approximate than the direct WKB approach of Salpeter 
and Van Horn and Schramm and Koonin (particularly, it neglects the curvature corrections). Numerically, 
Ichimaru, Ogata and lyetomi give a reaction rate which is close to the relaxed lattice model (model 2) of Salpeter 
and Van Horn The main reason for the coincidence of these rates is that the screening potential in the radial 
equation of Ichimaru, Ogata and lyetomi is close to the relaxed-lattice screening potential of Salpeter and Van Horn 
at ion separations r ~ 1.5 a, most important for pycnonuclear tunneling problem (see Fig. 2 in Ref. [s^). 

In spite of the differences in theoretical models 1-9, they result in similar reaction rates (Fig.lSJ. According to the 
above discussion, models 1 and 4 seem to be the most reliable among all available models for the bcc lattice, while 
models 5 and 8 seem to be the most reliable for fee. These reaction rates may be modified, for instance, by taking into 
account the quantum effect of the spreading of WKB trajectories or by a more careful treatment of the center-of-mass 
motion of reacting nuclei. Such effects will possibly reduce the reaction rate (as discussed in Ref. '53| with regard to 
the spreading of WKB trajectories). This could have been studied by direct PIMC simulations (e.g., Refs. _45^ 58j). 
PIMC is also a good tool to confirm the conclusions on dynamical effects of lattice response. However, PIMC is time 
consuming and requires very powerful computers. It is not clear whether today's computer capabilities are sufficient 
to obtain accurate PIMC pycnonuclear reaction rates. 

We suggest to calculate the reaction rates from Eq. H25|) taking into account that the constants Ccxp, Cp\ and Cpyc 
are not known very precisely. In particular, we propose two "limiting" purely phenomenologieal sets of these constants 
labeled as models 10 and 11 in Table These limiting parameters define the maximum and minimum reaction rates 
which enclose all model reaction rates 1-4 and 9 (proposed in the literature for the bcc lattice in a density range 
where the pycnonuclear carbon burning is important). They also enclose the most reliable models 5 and 8 for the fee 
lattice. 

The crucial parameter for modeling pycnonuclear fusion is the exponent Fpyc — exp(— Cpyc/VA) in Eq. I|23|l that 
characterizes the probability of Coulomb tunneling. It is easy to show that the exponent argument behaves as 
C'pyc/VA = a (ri2/rq„i)^ c>c p^^^®, where ri2 is the equilibrium distance between the interacting nuclei in their lattice 
sites, Tqm is the rms displacement of the nucleus due to zero-point vibrations in its lattice site, and a ~ 1 is a numerical 
factor which depends on a model of Coulomb tunneling. The usual condition is r^m ^ ''12 (and the tunneling length 
^ ''qm)- The exponent argument Cpyc/'v/A is typically large but decreases with growing p, making the Coulomb 
barrier more transparent. The tunneling is actually possible for closest neighbors (smallest ri2); the tunneling of 
more distant nuclei (higher ri2) is exponentially suppressed. Elastic lattice properties specify rqm and a, and are, 
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thus, most important for the reaction rate. The presence of different ion species, lattice impurities and imperfections 
may drastically affect the rate [l5l| . 

E. Thermally enhanced pycnonuclear regime 

The thermally enhanced pycnonuclear burning occurs with increasing T; it operates [l5j | in a relatively narrow 
temperature interval 0.5 Tp/ ln(l/VA) < 0.5 Tp. In this interval the majority of the nuclei occupy their ground 
states in a strongly coupled quantum Coulomb system, but the main contribution to the reaction rate comes from a 
tiny fraction of nuclei which occupy excited bound energy states. The increase of the excitation energy increases the 
penetrability of the Coulomb barrier, and makes the excited states more efficient than the ground state. 

The thermally enhanced pycnonuclear regime has been studied less accurately than the zero-temperature pycnonu- 
clear regime. Salpeter and Van Horn 'l5l calculated the thermally enhanced pycnonuclear reaction rate for models 
1 and 2 of a bcc lattice in the WKB approximation. The spectrum of excited quantum states was determined for a 
relative motion of interacting nuclei in an anisotropic harmonic oscillator field; the summation over discrete quan- 
tum states in the expression for the reaction rate was replaced by the integration. According to their Eq. (45), the 
enhancement of the reaction rate is approximately described by 

i?pyc(0) AV2 7^ VA )' ^ ^ 

where fi, fli, A, and Ai are model-dependent dimensionless constants. The exponent exp(— ATp/T) reflects the 
Boltzmann probability to occupy excited quantum states while the double exponent exp{(rii/\/A) e^^^^f describes 
the enhancement itself. In this case the characteristic energy of the reacting nuclei is 

Ep^^C,hijp + C2^ exp(-A,^y (27) 

where Ci and C2 are new dimensionless constants (^ 1). When T increases from T = to T ~ 0.5 Tp, the characteristic 
energy Epk increases from the ground state level, Epk ~ fuvp, to the top of the Coulomb potential well, -Epk ~ Z^e'^/a. 

The thermally enhanced pycnonuclear reaction rate was studied also by Kitamura and Ichimaru [53 adopting the 
formalism of Ogata, lyetomi and Ichimaru |5^ ("Sections IIII CI and IIII D|) . The relative motion of interacting nuclei 
was described by a model radial Schrodinger equation which employed the angle-averaged static MC plasma screening 
potential. The excited energy states were determined from the solution of this equation. Such an approach seems to 
be oversimplified. It gives the temperature dependence of the reaction rate (Eqs. (14) and (15) in Ref. [s^) which, 
functionally, differs from the temperature dependence, Eq. (|27|l . predicted by Salpeter and Van Horn. Nevertheless, 
numerically, both temperature dependencies at T < 0.5 Tp are in a reasonable qualitative agreement. 

We expect that the reaction rate in the thermally enhanced pycnonuclear regime will be further elaborated in the 
future. 

F. The intermediate thermo-pycnonuclear regime 

The intermediate thermo-pycnonuclear regime is realized at temperatures T ~ Tp (roughly, at Tp/2 < T < Tp) 
which separate the domains of quantum and classical ion systems. The main contribution to the reaction rate stems 
then from nuclei which are either slightly bound, or slightly unbound, with respect to their potential wells. The 
calculation of the reaction rate in this regime is complicated. We will describe this rate by a phenomenological 
expression presented in the following section. 

G. Single analytical approximation in all regimes 

Let us propose a phenomenological expression for the reaction rate which combines all the five burning regimes: 
R = Rpy,{T = 0) + Ai?(r), Ai?(T) = ^ S{Epi,) —^PF, 

T = exp(-f+C.f,e--./--A|), P = {fj ■ ^ 
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In this case, (/s = ^/V + F^]^/^; i?pyc(T = 0) is the zero-temperature pycnonuclear reaction rate (Section 
IIIID(I : AR{T) is the temperature-dependent part (with a product of an exponential function F and a pre-exponent 
P). The quantities r and T are similar to the familiar quantities t and T, but contain a "renormalized" temperature 
f: 



^ = ^2)' [f) ' ^-"^F' + (29) 

where Ct is a dimensionless renormalization parameter (^ 1). For high temperatures T ^ Tp we have t — > r, 

F — > F, and T ^ T. In this case the temperature dependent term tends to AR{T) — > Rth{T) ^ -Rpyc and Eq. H28(l 
reproduces the thermonuclear reaction rate fSections IIII Bl and IIII C|) . At low temperatures T <Tp the quantities r, 

F and T, roughly speaking, contain "the quantum temperature" Tp rather than the real temperature T in the original 
quantities r, F and T. In the limit of T we obtain F = 1/[VA (727r)i/6 Ct] and 7 = (s / (^^^^ C't ^ ) • 

At this point, let us require that at T ^ Tp the factor exp(— r) in the exponential function F, Eq. (|28|l . reduces to 
exp(— Coxp/\/A) in the exponential function Fpyc, Eq. (|23|l . This would allow us to obey Eq. (|26|l by satisfying the 
equality 

3^/¥/(2^/SC^/') = Cexp. (30) 

Taking Cexp we can determine Ct (see Table CJl. The double exponent factor in F, Eq. H28(l . will correspond to the 
double exponent factor in Eq. H2t)|) . Strictly speaking, Eq. 12t)|) contains two different constants A and Ai. However, 
taking into account the uncertainties of R in the thermally enhanced pycnonuclear regime (Section IIII E|) . we replace 
two constants by one. 

Finally, the quantity 7 in Eq. H28I) should be taken in such a way as to reproduce the correct limit 71 = 2/3 at 
T > Tp f Sect. IIII B|) and 72 = (2/3) (Cpi + 0.5) at T < Tp (see Eq. ^^). The natural interpolation expression for 7 
would be 

7= (T^i-f Tp272)/(T2 + Tp2). (31) 

In addition, we need the reaction energy iJpk to evaluate the astrophysical factor S'(£'pk). Since the S'-factor is a 
slowly varying function of energy, it is reasonable to approximate i?pk by the expression 

/ZV Tr\ / AT„\ , , 

Ep^^hu;p+i^-— + —j cxpi^—^j (32) 

which combines the expressions in the thermonuclear and pycnonuclear regimes. To avoid the introduction of many 
fit parameters (unnecessary at the present state of investigation), we set Ci = C2 = 1 in Eq. (|27|l . 
Thus, we propose to adopt the analytic expression (|28|l using the following parameters: 

(1) Csc = 1.0754 for the case of strong plasma screening in the thermonuclear regime fSection |lII C|l : 

(2) Cexp, Cpyc, and Cpi for conditions of zero-temperature pycnonuclear burning (see Tableland Section IlII D(l : 

(3) The quantum-temperature constant Ct fSection lIIID|l . which is important at T ^ Tp and expressed through Cexp 
via Eq. H30|) : the corresponding values of Ct are listed in Tabled 

(4) The last constant A, that is important at T '--^ Tp, is still free. 

We have checked (Fig. O that taking the optional model 1 from Table |2 and the value A = 0.5 results in a good 
agreement with the carbon burning rate calculated at p ~ 10^ — 10^° g cm^"^ and T < 0.5 Tp from Eq. (45) of Salpeter 
and Van Horn [l5l | (for the thermally enhanced pycnonuclear regime). (Notice that model 2 requires slightly lower 
A ^ 0.45.) Taking A = 0.35 leads to a noticeably higher reaction rate at T < 0.5 Tp, while taking A = 0.65 leads to a 
noticeably lower rate. 

Accordingly, for any model of zero-temperature pycnonuclear burning from TableQlwe suggest to adopt A — 0.5 as 
optional, A — 0.35 to maximize and A = 0.65 to minimize the reaction rate. In particular, model 1 with A = 0.5 seems 
to be the "most optional" ; our limiting model 10 from TableQlwith A = 0.35 is expected to give the upper theoretical 
limit for the reaction rate, while the other limiting model 11 with A — 0.65 is expected to give the lower theoretical 
limit. We also need the astrophysical factor S{E), which was described in Section ^1 for the carbon burning. We 
could easily introduce additional constants to tune our phcnomenological model when precise calculations of reaction 
rates appear in the future. 

For illustration. Fig. shows the temperature dependence of the carbon burning rate at p = 5 x 10^ g cm~'^. The 
solid curve is the most optimal model (based on both - zero-temperature and thermally enhanced - pycnonuclear 
burning models of Salpeter and Van Horn |T3| for the bcc static lattice). The double hatched region shows assumed 
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uncertainties of this model associated with variations of A from 0.35 to 0.65 (as if we accept the zero-temperature 
model but question the less elaborated model of thermal enhancement). The singly hatched region indicates overall 
uncertainties (limited by the models of the maximum and minimum reaction rates). The lower long-dashed line is 
obtained assuming classical thermonuclear burning without any screening fSection llll Bl) . The upper long-dashed line 
is calculated using the formalism of thermonuclear burning with screening fSection llll C|l . The screening enhancement 
of the reaction rate becomes stronger with the decrease of T . The formalism for describing this enhancement is 
expected to be valid at T > Tp, but we intentionally extend the upper long-dashed curve to T = 0.5 Tp, where the 
formalism breaks down and the curve diverges from the expected (solid) curve. The short-dash curve is calculated 
from the equations of Salpeter and Van Horn derived in the thermally enhanced pycnonuclear regime (model 1) 
and valid at T < 0.5 Tp. We intentionally extend the curve to higher T, where the formalism of thermally enhanced 
pycnonuclear burning becomes invalid and the curve diverges from the expected curve. Our phenomenological solid 
curve provides a natural interpolation at T ~ Tp between the short-dashed curve and the upper long-dashed curve. 

More complicated expressions for the reaction rate R in wide ranges of p and T were proposed by Kitamura [60j 
who combined the results of recent calculations of R in the different regimes. His expressions are mainly based on 
the results of Refs. which are not free of approximations (as discussed in Sections IHI CI IIII Dl and 

ini Ep . In contrast to our formula, Kitamura took into account the effects of electron screening (finite polarizability 
of the electron gas) and considered the case of equal and non-equal reacting nuclei. However, the electron screening 
effects are relatively weak; their strict inclusion in the pycnonuclear regime represents a complicated problem. We 
do not include them but, instead, take into account theoretical uncertainties of the reaction rates without electron 
screening. We have checked that the results by Kitamura (60l | for carbon burning in the most important T — p domain 
lie well within these uncertainties. 

Our formula gives a smooth behavior of the reaction rate as a function of temperature and density, without any 
jump at the melting temperature T = T^. We do not expect any strong jump of such kind since the liquid-solid 
phase transition in dense stellar matter is tiny. A careful analysis shows the absence of noticeable jumps of transport 
coefficients [6ll | and the neutrino emissivity owing to electron-nucleus bremsstrahlung. A direct example is given 
by the theory of nuclear burning. Ichimaru and Kitamura predicted a noticeable jump of the reaction rate at 
T = Tm, while a more careful analysis of Kitamura [6^] considerably reduced this jump. 

IV. CARBON BURNING AND IGNITION IN DENSE MATTER 

In this section we will analyze the rate of the ^^C-|-^^C reaction as a function of T and p and investigate the 
conditions for carbon burning in dense stellar matter. 

Because the probability for Coulomb tunneling depends exponentially on plasma parameters, changes in density p 
and temperature T have dramatic effects on the burning rate R. In thermonuclear regimes fSections IIII Bl and IIII C|) 
the ^^C-|-"'^^C rate is more sensitive to changes in temperature T than in density p. On the contrary, in pycnonuclear 
regimes (Sections IIII Dl and IIII E|) the rate depends significantly on the density p. For instance, if T decreases from 
3 X 10^ K to 3 X 10^ K at p = 5 X 10^ g cm~^ (Fig. El thermonuclear burning with strong screening), the reaction 
rate drops by ^ 20 orders of magnitude. Neglecting the enhancement due to plasma screening, the rate will drop 
by ten more orders of magnitude. An increase in density p from 10^ g cm~'^ to 10^"'^ g cm~'^ at T < 3 x 10^ K (in 
the zero-temperature pycnonuclear regime) results in a rate increase of ~ 100 orders of magnitude (Fig. |3J). Note 
that no carbon can survive in a degenerate matter at p > 3.90 x 10^° g cm~^ because of the double electron capture 
^'^C ^|B "'^fBe (e.g., Shapiro and Teukolsky H^). The electron capture has a well defined density threshold, 
3.9 x 10^° g cm~^, and proceeds quickly after the threshold is exceeded. We will ignore this process in the present 
section. 

The strong dependence of the rate R on density p and temperature T leads to huge variations of the characteristic 
time scale Tbum = rii/R for carbon burning. Figure^lshows two solid lines in the p — T plane, along which Tburn = 1 s 
and lO^'' years (nearly the Universe age), respectively. They are calculated using the most optional carbon burning 
model (model 1 from Tabled A — 0.5). The lines are almost horizontal in the thermonuclear burning regime {R is 
a slowly varying function of p) and almost vertical in the pycnonuclear regime (i? is a slowly varying function of T) . 
The bending part of the lines corresponds to the thermally enhanced pycnonuclear and intermediate thermo-pycno 
nuclear regimes. At T and p above the upper line the burning time is even shorter than 1 s; at these conditions no 
carbon will survive in the dense matter of astrophysical objects. For conditions below the lower solid line, Tbum is 
longer than 10^° years, and carbon burning can be disregarded for most applications. 

Thus, the studies of carbon burning can be focused on the narrow strip in the p — T plane between the lines of 
Tburn = 1 s and Tburn = 10^" yr. Hatched regions show theoretical uncertainties of each line (limited by the maximum 
and minimum reaction rate models, Section flll Cp . The uncertainties are relatively small in the thermonuclear regime 
where all models give nearly the same reaction rate. The uncertainties are higher in other burning regimes. 
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Having a carbon burning model we can plot the carbon ignition curve. This curve is a necessary ingredient for 
modeling nuclear explosions of massive white dwarfs (producing supernova la events, so important for cosmology; 
see, e.g., Refs. 0, |^) and for modeling carbon explosions of matter in accre ting neutron stars (viable models of 
superbursts observed recently from some accreting neutron stars; e.g., Refs. 

The ignition curve is commonly determined as the line in the p—T plane (Fig.^ where the nuclear energy generation 
rate is equal to the local neutrino emissivity of dense matter (the neutrino emission carries the generated energy out 
of the star). At higher p and T (above the curve) the nuclear energy generation rate exceeds the neutrino losses and 
carbon ignites. In Fig. ^ we present the carbon ignition (solid) curve, calculated using the most optional model of 
carbon burning, together with its uncertainties (limited by the minimum and maximum rate models). The neutrino 
energy losses are assumed to be produced by plasmon decay and by electron-nucleus bremsstrahlung. The neutrino 
emissivity due to plasmon decay is obtained from extended tables calculated by M. E. Gusakov (unpublished); they 
are in good agreement with the results by Itoh et al. [65l |. The neutrino bremsstrahlung emissivity is calculated using 
the formalism of Kaminker et al. j62|, which takes into account electron band structure effects in crystalline matter. 

For p < 10^ g cm^^ theoretical uncertainties of the ignition curve are seen to be small. They become important at 
P ^ 10^ g cm~^ and T ~ (1 — 3) x 10® K in the intermediate thermo-pycnonuclear burning regime and the thermally 
enhanced burning regime. This p — T range is appropriate for central regions of massive and warm white dwarfs which 
may produce type la supernova explosions. Lower T are also interesting for these studies (e.g., Baraffe et al. 0)- 

If we formally continue the ignition curve to lower T, it will bend and shift to lower densities, where the nuclear 
burning time scale Tburn is exceptionally slow exceeding the age of the Universe. The bend is associated with a very 
weak neutrino emission at T < 10® K. These parts of the ignition curve are oversimplified because at low T the 
energy outflow produced by thermal conduction becomes more efficient than the outflow due to the neutrino emission. 
These parts are shown by the long-dash line (and their uncertainties are indicated by thin dash-and-dot lines). 
Unfortunately, the conduction energy outflow is non-local and "non-universal" . It depends on specific conditions 
of the burning environment (a white dwarf core or a neutron star crust) and the associated thermal conductivity 
(provided mainly by strongly degenerate electrons). In this case the ignition becomes especially complicated. A very 
crude estimate shows that the ignition curve, governed by the thermal conduction, is nearly vertical and close to the 
Tburn = 10^" yr curve in the range of T from 10® K to 10^ K in Fig. El At T < 10^ K the curve is strongly affected 
by the thermal conductivity model. In a cold ideal carbon crystal, umklapp processes of electron-phonon scattering 
are frozen out (e.g., Ref. 66]). Under these conditions the electron conduction is determined by inefficient normal 
electron-phonon scattering, leading to high conductivity values. This shifts the ignition to higher p. On the other 
hand, carbon matter may contain randomly located ions of other elements (charged impurities) which can keep the 
electron Coulomb scattering rather efficient and maintain a low electron thermal conductivity. In this case the ignition 
curve at T < 10^ K remains nearly vertical. 

V. CONCLUSION 

The goal of this paper was to develop a phenomenological formalism for calculating fusion reaction rates between 
identical nuclei. This formalism should be applicable for a broad range of thermonuclear and pycnonuclear burning 
scenarios. It involves a generalized treatment for calculation of the fusion probability at low energies and the devel- 
opment of a single simple phenomenological expression for the fusion rate valid in a wide range of temperatures and 
densities. 

We have introduced a generalized model approach for calculating the S'-factor of heavy-ion fusion reactions relevant 
for stellar nucleosynthesis processes. We have demonstrated the applicability and reliability of the approach by 
calculating the astrophysical factor S{E) for the carbon fusion reaction ^^C+^^C (Sectionnj and by comparing the 
theoretical results with experimental data. 

Furthermore, we have analyzed fSection IIII|) previous calculations of the fusion rate for identical nuclei in stellar 
matter, with emphasis on the complicated problem of Coulomb barrier penetration in a dense-plasma environment. 
Combining the results of previous studies, we have proposed a single simple phenomenological expression for the 
fusion rate, valid in all five fusion regimes (that can be realized in the different p—T regions). Our formula contains 
adjustable parameters whose variations reflect theoretical uncertainties of the reaction rates. 

For illustration, we have considered fSection II VI) the efficiency of carbon burning in dense matter and the conditions 
for carbon ignition in white dwarf cores and neutron star crusts. We show that carbon burning is actually important 
in a sufficiently narrow p — T strip which is mainly determined by the temperature T ~ (4 — 15) x 10® K as long 
as yO < 3 X 10^ g cm^^, and by the density p ~ (3 — 50) x 10^ g cm~'^ as long as T < 10® K. On the basis of these 
results we suggest that the current knowledge of nuclear fusion is sufficient to understand the main features of carbon 
burning in stellar matter, especially at p < 3 x 10^ g cm~^. 

We have focused on the simplest case of heavy-ion burning in a one-component Coulomb system; particularly, in 
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a perfect crystal. There is no doubt that dense matter of white dwarfs and neutron stars are more complicated and 
require a more complex approach taking into account mixtures of different heavy nuclei and imperfections in dense 
matter. The complexity ranges from essentially two-component plasma conditions anticipated in the carbon-oxygen 
cores of white dwarfs to the multi-component isotope distribution in the ashes of accreting neutron stars [g^ . 

In a forthcoming paper we will expand the presented analysis to the case of the fusion rates between different 
isotopes. We will employ this formalism for calculating the iS-factors for a broad range of heavy-ion fusion reactions. 
We will include the results in a pycno-thermonuclear reaction network and simulate the nucleosynthesis in high density 
stellar matter. 
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TABLE L The table presents the coefBcients Coxp, Cpyc and Cpi of the pycnonuclear reaction rate obtained at T = (see 
I25|l . The dimensionless parameter Ct, which is related to the "renormalized" temperature (Eq. 129^ . is also included in 
table. 



No. 




Cpyc 


Cpi 


Ct 


Model 


Refs. 


1 


2.638 


3.90 


1.25 


0.724 


bcc; static lattice 




2 


2.516 


4.76 


1.25 


0.834 


bcc; relaxed lattice - WS 




3 


2.517 


4.58"' 


1.25 


0.834 


bcc; relaxed lattice 


\16} 


A 




U. lO 


1 9^ 




Dec, eiicCLiVc maab appiox. 




5 


2.401 


7.43"' 


1.25 


0.960 


fee; static lattice 


m 


6 


2.265 


13.5"' 


1.25 


1.144 


fee; relaxed lattice - WS 


m 


7 


2.260 


12.6"' 


1.25 


1.151 


fee; relaxed lattice 


m 


8 


2.407 


13.7"' 


1.25 


0.953 


fee; effective mass approx. 


fl6] 


9 


2.460 


0.00181 


1.809 


0.893 


bcc; MC calculations 


[53] 


10 


2.450 


50 


1.25 


0.904 


maximum rate 


present paper 


11 


2.650 


0.5 


1.25 


0.711 


minimum rate 


present paper 



Corrected for the curvature factor as explained in the text. 
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FIG. 2: Astrophysical factor S{E) as a function of the center-of-mass energy E, derived from experimentally measured 
cross sections. Lines show theoretical results obtained within the barrier penetration model for the different model density 
distributions (see text for details). Various symbols present experimental results. 
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FIG. 3: Elastic scattering angular distributions for the C+ C system at energies around and slightly above the Coulomb 
barrier jlj . The lines are the results of an optical model calculation assuming the Sao Paulo potential to describe the real part 
of the nuclear interaction, combined with a simple model to describe the imaginary part of the optical potential (see details in 
the text). 
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FIG. 4: Temperature-density diagram for carbon matter. Short-dashed lines show the electron degeneracy temperature Tp, 
the temperature T; of the appearance of ion liquid, the melting temperature Tm of ion crystal, and the ion plasma temperature 
Tp. Solid lines correspond to the carbon burning times Tbum = 1 s and rbum = 10^'' yr, and to carbon ignition; they are 
calculated using the most reliable model of carbon burning (Section IIII Gl . Hatched strips show theoretical uncertainties of 
these lines (limited by the minimum and maximum reaction rate models). The long-dashed line exhibits the unreliable part of 
the ignition curve; nearby thin dashed-and-dot lines (to the right and left) indicate its assumed uncertainties. 
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FIG. 5: Rate of pycnonuclear carbon burning at T = as a function of density for the different theoretical models (from Table 
01. Solid and dashed lines refer to the burning in bcc and fee crystals, respectively; the dash-and-dot line 'K' is the model by 
Kitamura |60ll (for bcc crystal). Hatched strip shows assumed uncertainties of the reaction rates for bcc crystals (limited by 
models 10 and 11 from Table 0. 



22 




- pycnonuclear /thermonuclear- 

I I LJ I I I I I I I / I I Lj LJ I I I I I l_ 

7 7.5 8 8.5 9 9.5 



logio T [K] 



FIG. 6: Temperature dependence of the carbon fusion rate at p = 5 x 10^ g cm"''. The sohd hne 1 is our most optional 
interpolation expression (Sect. HTl GL based on model 1 from Table |I] with A — 0.5. Doubly hatched region shows theoretical 
uncertainties of model 1 associated with variations of A from 0.35 to 0.65. The dash-and-dot line 'K' is the interpolation of 
Kitamura [g^- The short-dashed line 1 is calculated from the expressions of Salpeter and Van Horn Hsf . which are valid in 
the pycnonuclear regime (T = and with the thermal enhancement). Long-dashed lines show the thermonuclear reaction rates 
calculated with account for plasma screening (Sect. ITTl Cll and without screening (Sect. HTlBII . Singly hatched region displays 
total assumed theoretical uncertainties of the reaction rates. 



